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Abstract 


A two-dimensional steady state model for a PEM fuel cell cathode is described in this work. All the components in the cathode such as the 
gas manifold, diffusion layer, microporous layer and the catalyst layer are modeled. The effect of the liquid water is taken into account in every 
layer of the cathode. The model was developed and simulated using a combination of Maple and MATLAB. The combination provides a flexible 
framework for quickly developing models with various assumptions and different complexities. The cathode catalyst layer was modeled using 
both macrohomogeneous and spherical agglomerate characterizations. The model is validated using experimental data. During model validation, 
various assumptions are considered for transport within the porous layers of the cathode. Subsequently, the assumptions and characteristics that 
best predicts the experimental data are highlighted. The major conclusion of this work is that a model that includes liquid water in all the layers with 
a flooded spherical agglomerate characterization for the reaction layer best predicts the PEM fuel cell behavior in terms of an i—v characterization 
for a wide range of reactant flow rates. The utility of the steady state model for the optimization of the cathode catalyst layer design parameters is 


also described. 
© 2007 Published by Elsevier B.V. 
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1. Introduction 


There is considerable interest in the modeling of proton 
exchange membrane fuel cells (PEMFC), and a number of 
PEMFC models have been proposed in the literature over the 
last two decades. Models proposed during the early years were 
typically one-dimensional and accounted for steady state mass 
transport and electrochemical kinetics. Subsequently, both sim- 
plified and complex models in terms of dimensionality and 
physicochemical phenomena have been studied. These models 
have been used for a variety of purposes such as, prediction of the 
typical characteristic (current—potential) curves, parameter and 
operating conditions sensitivity analysis, and three-dimensional 
temperature, pressure and species concentration distributions in 
the case of fuel cell stacks. 
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There are two main approaches that have been pursued in the 
modeling of PEMFCs. One is the detailed models of transport 
processes and electrochemical reactions that take place in a fuel 
cell. The other approach is the development of simplified models 
with well defined reactor conditions for correlation of fuel cell 
operation as exemplified by the work of [1]. In this paper, we 
follow the former approach. Our aim here is the development of 
a model that can be used for detailed analysis and optimization 
of fuel cell systems. This would entail including the effect of 
liquid water in all the layers of the fuel cell. Further, we also 
require that the experimental parameters to be correlated to the 
fuel cell performance. Development of such a detailed model- 
ing and computational approach would allow for comprehensive 
phenomenological study of several important factors as outlined 
excellently in [2]. 

A review of state-of-the-art in models proposed for PEM 
fuel cells is presented in the next section. A detailed review of 
the present status of fundamental models for fuel cell engineer- 
ing is also presented in [3]. The review highlighted the current 
status of hydrogen/air polymer electrolyte fuel cells (PEFCs), 
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S ls L length of the gas flow channels (m) 
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Dieffm effective diffusivity of species i in the microp- of the catalyst layer (mol m~? s~') 
orous layer (m? s7!) Rw interfacial transfer of water between liquid and 


vapor (mol m~? s7!) 
RHair relative humidity of cathode inlet air (%) 
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Tair cathode inlet air temperature (K) 
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a transfer coefficient 
ômem thickness of ionomer film covering the agglomer- 
ate (m) 
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Ow thickness of water layer on top of the agglomerate 
(m) 

Ed void fraction inside the gas diffusion layer 

Em void fraction inside the microporous layer 

€k void fraction inside region k 

€mem fraction of volume occupied by the ionomer inside 
the catalyst layer 

Er void fraction inside the catalyst layer 
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Nr local overpotential inside the catalyst layer (V) 

é contact angle 

Keff effective conductivity of ionomer inside the cata- 
lyst layer (mho m7!) 

Kmem Conductivity of ionomer (mho m7!) 

Os oxygen excess ratio 

Àw water content inside the ionomer (mol 
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Uw viscosity of liquid water (kg m7! s7!) 

Pe density of carbon (kg m~?) 
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pr local ionomer potential inside the catalyst layer 
(V) 

Yy Thiele modulus 

Subscripts 

i index for the species: O2, N2, H2O 

k index for the region: gas flow channel, diffusion 
layer, microporous layer and catalyst layer 


direct methanol fuel cells (DMFCs), and solid oxide fuel cells 
(SOFCs). Our main observation on the detailed models available 
in the literature is that most of these models do not characterize 
the effect of liquid water on the fuel cell performance. Com- 
prehensive models that include the effect of liquid water in our 
view are the ones proposed by [4-6]. Pasaogullari and Wang 
[4] describe the governing physics of water transport in both 
hydrophilic and hydrophobic diffusion media along with one- 
dimensional analytical solutions of related transport processes. 


Higher-magnification SEM image, showing 


much smaller particles 


Still higher-magnification SEM image, 
that the particles are in fact agglomerates of Showing that the particles are agglomerates 


of much smaller particles 


They report that the liquid water transport across the gas dif- 
fusion layer (GDL) is controlled by capillary forces resulting 
from the gradient in phase saturation. A one-dimensional ana- 
lytical solution of liquid water transport across the GDL was 
derived. Effect of GDL wettability on liquid water transport was 
explored in detail for the first time. Furthermore, the authors 
also investigate the effect of flooding on oxygen transport and 
cell performance and show that flooding diminishes the cell per- 
formance as a result of decreased oxygen transport and surface 
coverage of active catalyst by liquid water. Lin et al. [5] model 
the liquid water in the gas diffusion and catalyst layers and they 
characterize the catalyst layer through a cylindrical geometry. 
Their model domain consists of the membrane, the cathode cat- 
alyst layer and the cathode diffusion layer. Further, experimental 
studies have shown that the catalyst layer in PEMFC is better 
characterized by spherical agglomerates, see Fig. 1[7]. In the 
work of [6], liquid water is modeled in all the layers including 
the manifold. However, the catalyst layer is not characterized 
and a macrohomogeneous approach is used. We later show with 
experimental data that this would lead to a poor prediction of 
fuel cell behavior in certain operational regimes. 

In this work, the first major contribution is a model that 
includes liquid water in all the layers of the cathode (gas man- 
ifold, gas diffusion layer, microporous layer and the catalyst 
layer) and the catalyst layer is characterized as spherical agglom- 
erates, in line with the experimental evidence [7—9]. Further, we 
compare several models with different simplifications and show 
the importance of modeling the liquid water in all the layers with 
the correct characterization of the catalyst layer. Cathode exper- 
imental data for a wide range of flow rates is used to validate and 
evaluate the proposed models. Our main conclusion is that when 
the catalyst layer is characterized using spherical agglomerates 
with liquid water effects, and the transport in the porous regions 
is by diffusion in the bulk, the model predictions match very well 
with experimental data for the entire polarization range and at 
various flow rates. 

The second major contribution of this work is the correlation 
of the experimental variables used in the preparation of the 
membrane electrode assembly (MEA) to the model parameters. 
This correlation is derived based on a spherical agglomerate 
characterization of the reaction layer. Hence, once the exper- 
imental variables used in the preparation of the MEA such as 


Higher electron energy in the 
HR-SEM shows that the agglomerates 
are built up from 30 nm carbon 
particles, and that these particles are 
coated with finely dispersed platinum 


Fig. 1. Scanning electron microscope (SEM) images of PEM fuel cell electrode [7]. 
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the thickness of the catalyst layer (trL), the catalyst loading 
(mp), the weight fraction of platinum on carbon (fpr), the 
weight fraction of ionomer inside the catalyst layer (fmem), 
and the density of ionomer (Pmem) are provided, all the other 
model parameters are derived through balance equations. For 
example, the void fractions and the thickness of the ionomer 
layer covering the agglomerates are related to the experimental 
parameters. This directly ties in the model predictions to the 
experimental characterization. We believe that this is of major 
importance because directly usable multi-parameter optimiza- 
tion studies are possible with this model. The results of such an 
optimization study have been reported in another paper [10]. 


2. Literature review 


Recent literature is abundant with a variety of PEMFC 
models. Both simplified and complex models in terms of dimen- 
sionality and physicochemical phenomena have been studied 
[11-43]. Seminal papers in PEMFC modeling were published 
by Springer et al. [36] and Bernardi and Verbrugge [11,12] and 
in contemporary literature, significant contributions were made 
by [4,44]. A common feature in most of these models is that the 
reaction or catalyst layer is not modeled in detail. The reaction 
layer is treated as an ultra-thin layer, thus neglecting the trans- 
port of reactant gases and products. Hence, the catalyst layer is 
treated as a source/sink boundary condition for transport equa- 
tions in the gas diffusion layer. Contrary to this assumption, 
even if gas phase transport is neglected on the consideration of 
an ultra-thin layer, the presence of ionomer in the reaction layer 
along with carbon and platinum makes transport within the pores 
of the ionomer important. Moreover, catalyst layer is the region 
where various limiting mechanisms can occur and thus, can have 
a strong influence on the overall performance of the cell. 

As gas diffusion electrodes (GDEs) are difficult to character- 
ize, one of the first assumptions that was made to model them 
was the concept of “flooded agglomerates”, introduced by Giner 
and Hunter [45]. They have considered cylindrical geometry for 
the agglomerates. Results were presented for alkaline oxygen 
electrode. The potential drop was assumed to change only in the 
axial direction and diffusion of the dissolved gases was assumed 
to be in the radial direction of the agglomerates. The effect of the 
cylindrical agglomerate radius on the current generated and its 
distribution were studied. Porous GDEs with the same assump- 
tions have been extended to model the phosphoric acid fuel cells 
(PAFC) cathode and anode in detail [47,46]. Various transport 
and kinetic processes which take place in the porous electrodes 
were taken into account. The model was used in the simula- 
tion mode for predictive analysis and for electrode development 
process. 

One of the drawbacks with the above-proposed flooded- 
agglomerate model is that it does not consider any tortuosity 
for the gas phase transport as the agglomerates are assumed 
to completely extend from the gas side to the electrolyte side. 
The cylindrical flooded-agglomerate model was modified by 
Celiker et al. [48] considering spherical geometry. They have 
investigated their model predictions by considering the cathodic 
reduction of oxygen in alkaline medium. Subsequently, many 


studies conducted by various researchers with the spherical 
flooded-agglomerate model were presented for alkaline fuel 
cells (AFC) [49,50] and PAFC [50-52]. 

Even in the case of PEM fuel cells, researchers have studied 
the effect of various phenomena in the catalyst layer based on 
flooded-agglomerate model [50,53-57]. Perry et. al. [50] have 
developed a model for gas diffusion electrode that can be used 
as a diagnostic tool for designing of fuel cells. This is a one- 
dimensional model for mass transport in the zone where Tafel 
kinetics is valid. The models presented were generally valid 
for any GDE with either liquid electrolyte (AFC and PAFC) 
or ion-exchange membrane (PEMFC). The model was used to 
study the effects of mass-transport limitations on the polariza- 
tion characteristics of oxygen reduction reaction in the cathode. 
Using negligible mass transfer resistance in the gas phase the 
model develops a function for evaluating the Thiele modulus 
for the catalyst-binder agglomerates. These relations along with 
ion transport equations are combined to develop a single variable 
second order differential equation. For this equation, asymptotic 
solutions were developed at different limiting conditions. The 
model also predicts different Tafel slopes for distinct regions. 
The authors have also shown how the results may be used as 
a diagnostic tool for analyzing fuel cell cathode data. Siegel et 
al. [54] have proposed a steady state two-dimensional PEMFC 
model based on agglomerate geometry for the catalyst layer. 
The agglomerates are characterized by mean diameter and a 
characteristic length. Based on the model results, it has been 
highlighted that the fuel cell performance is highly dependent 
on catalyst structure. Wang et al. [56] have investigated trans- 
port and reaction kinetics in spherical agglomerates of cathode 
catalyst layer. They have considered two types of spherical 
agglomerates: the first one consisting of a mixture of car- 
bon/catalyst particles and perfluorosulfonated ionomer (PFSI) 
and the other type consisting of carbon/catalyst particles and 
water-filled pores. The model has been used to study current 
conversion, reactant and current distribution and catalyst utiliza- 
tion. However, most of these models do not treat liquid water in 
the catalyst layer. Further these models ignore the other layers 
in an electrode. 

One of the probable reasons for neglecting reaction layer 
in PEMFC models published in the beginning could be lack 
of instrumentation to characterize its morphology accurately. 
With the availability of advanced microscopy instruments like 
scanning electron microscope (SEM) and transmission elec- 
tron microscope (TEM), researchers have been able to study 
the morphology of complex nanostructures such as, PEM 
fuel cell electrodes. Middleman [7] has studied the structure 
of membrane-electrode-assembly (MEA) using high-resolution 
scanning electron microscopy (HR-SEM). In his investigation 
he has shown that the catalyst layer consists of a random distri- 
bution of pores and particles, see Fig. 1. It was also shown, using 
higher magnification, that the particles are agglomerates of much 
smaller particles coated with a film of Nafion. From the images in 
the figure it can be clearly seen that the agglomerates are spher- 
ical in shape. Lee et al. [8] and Liu et al. [9] have also published 
their investigations of PEM fuel cell electrodes using SEM/TEM 
that corroborate Middleman’s work. Therefore, spherical 
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agglomerates can be treated as realistic representation of cathode 
catalyst layers in PEM fuel cells for modeling purposes. The util- 
ity of this conceptualization over macrohomogeneous character- 
ization is clearly described during model validation, Section 5. 


3. Steady state model 


A schematic of the PEM fuel cell cathode that is modeled in 
this work is illustrated in Fig. 2. The schematic shows the bipo- 
lar plate at the top, Toray carbon paper (TGP 120) as a diffusion 
layer, microporous layer below that, and the catalyst layer at the 
bottom. The bipolar plate consists of straight channels of uni- 
form cross-section. As air travels along the channels, it transports 
through the diffusion layer and microporous layer and enters the 
catalyst layer. Typically, fuel cell catalyst layer models are based 
on the assumption of either macrohomogeneous or flooded- 
agglomerate structure [21,45,49,51,54,58,59]. In this work, the 
catalyst layer has been modeled considering both the macroho- 
mogeneous and the spherical flooded-agglomerate structure. For 
the case of PEM fuel cells, a flooded-agglomerate is a uniform 
mixture consisting of Pt nanoparticles supported on carbon with 
the hydrated ionomer in the micropores. Hence, the catalyst layer 
consists of a cluster of flooded-agglomerates with free space in 
between for the gas to diffuse through and reach the surface of 
each flooded-agglomerate. In addition, each spherical agglom- 
erate is assumed to be coated with a thin film of ionomer. Fig. 3 
shows the schematic of the catalyst layer and an enlarged view of 
the spherical agglomerate. The figure also shows the membrane 
layer below the catalyst layer. At the surface of the ionomer, 
oxygen present in the gas dissolves into the water present inside 
the pores of the ionomer. At the surface it is assumed that there 
exists an equilibrium between the partial pressure of oxygen 
in the gas phase and the dissolved concentration in the ionomer 
phase. The dissolved oxygen diffuses through the ionomer pores 
and reaches the active catalyst sites, where the following oxygen 
reduction reaction takes place: 


O2 + 4H+ +4e7 > 2H,0 (1) 


Gas Diffusion 


Catalyst Layer 


Diffusion Layer 
Gas Inlet 


Fig. 2. Schematic of a PEM fuel cell cathode. 


The hydrogen ions produced in the anode catalyst layer travel 
across the membrane layer and reach active catalyst sites inside 
the cathode through a network of micropores in the ionomer. 
Detailed model equations for the gas flow channel, diffusion 
layer, microporous layer and catalyst layer of the PEM fuel cell 
are described below. The reader is referred to the nomencla- 
ture for details on variables and parameters description. The 
following assumptions are considered for setting up the model 
equations: 


Assumption 1. Isothermal conditions are considered through- 
out the region of interest. 


Assumption 2. Pressure gradients in the X direction in all the 
regions are negligible. Hence, velocities in the X direction are 
zero. 


Assumption 3. Effective diffusivities are assumed for diffusive 
transport in the gas phase inside the porous regions. 


Assumption 4. Butler—Volmer kinetics are considered for the 
oxygen reduction reaction. 


Assumption 5. Physical properties of the ionomer inside the 
catalyst layer are considered same as that of the membrane. 


Assumption 6. Potential drop in the solid phase due to resis- 
tance to the electron transport is assumed to be negligible. 


Assumption 7. The gas mixture inside the region of interest is 
assumed to behave as an ideal gas. 


Even though some of the above assumptions are generally 
not valid for all cases, they have been considered for various 
reasons. The assumption of isothermal conditions was based on 
the fact that the unit cell considered for model validation was 
small, with an active area of 20cm. For more detailed two- 
phase studies, it will be important to consider non-isothermal 
effects in the various regions of the fuel cell [39]. Meng and 
Wang [60] investigated the effects of electron transport through 
the gas diffusion layer (GDL) for the first time. They show that 
the current distribution was determined by two factors: oxygen 
supply and lateral electronic resistance in GDL. At a high cell 
voltage, the lateral electronic resistance dictated the current dis- 
tribution and at low cell voltages, oxygen concentration played a 
dominant role in determining the current distribution. However, 
potential drop effects due to electron transport were assumed 
negligible in order to avoid making the model more complex 
and not to deviate from the main focus of the paper—to high- 
light the liquid water effects in the different regions of the fuel 
cell cathode. 


3.1. Gas flow channel 


The following assumptions are considered for the gas flow 
channels: 


e All channels in the graphite plate are assumed to have equal 
air flow rate. 

e Negligible edge effects. This would result in uz = 0 

e Based on the previous assumption,u, = 0 


380 R.M. Rao et al. / Journal of Power Sources 173 (2007) 375-393 


Membrane Layer 


Hydrated 
Ionomer 


Platinum 
Nanoparticles 


Carbon 


Single agglomerate 


Fig. 3. Schematic of a PEM fuel cell catalyst layer and spherical agglomerate. 


e The boundary condition of no slip at the walls has also been 
relaxed and a constant velocity in Y direction is assumed. This 
gives us 


Uy = Uinlet (2) 


This simplification has been made in spite of the existence of 
a fully developed velocity profile for rectangular ducts [61] 

e The liquid water inside the gas flow channels is assumed to 
exist in the form of tiny droplets and travel with the gas veloc- 
ity [62]. Hence, mist flow model is applied to describe liquid 
water transport in the gas flow channels 


Hence, gas flow and liquid water transport in the channels 
can be considered as a plug flow with simultaneous exchange of 
species at the boundary between the gas flow channels and the 
backup substrate. A schematic of a control volume inside the 
gas flow channel is illustrated in Fig. 4. 

The species molar balance equation inside the gas flow chan- 
nels can be written as 


ð 
— — (Ci, gcUinlet) -V.-Ji=0 (3) 
dy 
Liquid water 
Y=0 droplets Control volume Y=L 
| 
Y 
= reo. Excess feed 
X A, + 
| eo 
Feed input Bulk flow “tog! > Product gas 
— io ! —> Naina S 


ofgas ~~ } oo: 


ist -t-4 


Exchange between 
gas channel and diffusion layer 


Fig. 4. Control volume inside gas flow channel. 


where, J; is the flux due to diffusion of species / into the diffusion 
layer from the gas flow channel. It is calculated from the local 
flux evaluated at the boundary x = 0 and can be expressed as 


Ji = — Di effa V Cialx=o (4) 


where, Dj effa is the effective diffusivity of species i inside the 
diffusion layer and C; q is the concentration of species i inside 
the diffusion layer. 

For water vapor, an additional term accounting for evapora- 
tion/condensation appears in the above species balance equation. 
Hence, the conservation equation for water vapor can be written 
as 


f) 
~ gy CHaO0.getinler) - V - Jmo — Ry =0 (5) 


where Rw is the interfacial transfer of water between liquid and 
water vapor and is defined as [24] 


e(l — s) 
R = k sat 
w c RTocl Yw(Pw Pw Jq 
EKSP 
+ ky——(pw — pN — q) (6) 
My 


where ke and ky are the condensation and evaporation rate con- 
stants, respectively; eg is the void fraction (= 1.0 for gas flow 
channel); s is the liquid water saturation level, which is the frac- 
tion of void volume occupied by liquid water; yy and py are 
mole fraction and partial pressure of water vapor, respectively; 
and p**' is the water vapor saturation pressure. The parameter q 


is a switching function that is defined by [24] 
j= 1.0 + (|Pw — Px )/(Pw = es 
E 2 


Hence, if pw > p$* then q = 1 and the interfacial transfer of 
water will be a condensation process. On the other hand, if pw < 


(7) 
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ps" then q = 0 and the interfacial transfer of water will be an 
evaporation process. 

The conservation for liquid water in the gas flow channel can 
be written as 


ð 
My By Cene -V. Nwa + Ry =0 (8) 


where Sg¢ is the fraction of volume occupied by the liquid water 
in the gas flow channel; Nw,q is the flux of liquid water into the 
gas flow channel from the diffusion layer. It is given by Darcy’s 
law and is described in the next section. 


3.2. Diffusion Layer 


Since no pressure gradients are assumed inside diffusion 
layer, transport of species is governed purely by diffusion due 
to concentration gradients. Moreover, in order to simplify com- 
putations, diffusive transport is defined using Fick’s law with 
an effective diffusivity instead of using Stefan—Maxwell equa- 
tions. Recent numerical studies [63,64] have shown that Fick’s 
formula with effective diffusion coefficient leads practically to 
the same results as those with Stefan—Maxwell equations. In the 
absence of any reaction in the diffusion layer, species conserva- 
tion equations can be written as 


=V + (—DyettaV Cia) — Rw = 0 (9) 


where, Dj effa is the effective diffusivity and C; q is the con- 
centration of species i inside the diffusion layer. The effective 
diffusivities of species in porous regions are related to the dif- 
fusivities in gaseous regions by the equation 


Diett = &! °C — s)? Dim (10) 


where eg is the porosity of the region k, s is the liquid water sat- 
uration and Dim is the diffusivity of the species i in the mixture, 
which is related to the binary diffusivities and is given by [65] 


-1 
N 


D an 


j=j 


Dim = (l — yi) 


The binary diffusion coefficients in the above equation are esti- 
mated using Chapman—Enskog formula [66]. The term Rw in 
Eq. (9) is applicable only for water vapor conservation, which 
defines the volumetric rate of evaporation/condensation and is 
given in Eq. (6). 

In the porous regions ( diffusion layer, microporous layer 
and catalyst layer), the liquid water is driven by capillary force. 
Darcy’s law is used to describe the flow of liquid water 


_ Ky(s) 


w 


u = V(Pi) (12) 


Therefore, the molar flux of liquid water can be written as 


PwKw(s) 


Ny =- 
v My bw 


V(P) (13) 


where py, My, and uw are density, molecular weight, and vis- 
cosity of liquid water, respectively; Ky,(s) is the permeability 


of liquid water in the porous regions. The capillary pressure is 
defined as the difference between total gas pressure (P,) and 
pressure of liquid water (P4) 


Po = Pg — Pi (14) 
Therefore, Eq. (13) can be written as 
PwKwy(s) 
Ny = V(Pe — P, 15 
w Melie (Pg — Pc) (15) 


Since gas phase pressure gradients are assumed negligible inside 
porous regions, the above equation reduces to 


_ PwKwy(s) 


Ny = Wis V(—Pe) (16) 
_ PwKw(s) dP. 
= Mali ( ae ) Vs (17) 


Both capillary pressure (P,) and permeability (Ky(s)) are func- 
tions of liquid water saturation (s). Several empirical expressions 
are available to describe the dependence of P, and Ky(s) on the 
liquid water saturation [24,62]. In order to reduce the number of 
fitting parameters in the model, (—dP,/ds) is treated as a con- 
stant and Ky(s) is assumed to depend linearly on liquid water 
saturation, i.e., Ky(s)= Kwos, where Kwo is the permeability of 
liquid water at 100% saturation [53]. 

Writing a conservation equation for liquid water in the 
absence of any reaction in the diffusion layer leads to 


-V -Nyw+ Ry =0 (18) 
3.3. Microporous layer 


For modeling purposes, the microporous layer is similar to the 
diffusion layer. The only difference is in the physical parameters 
such as, the void fraction, pore structure and thickness. Hence, 
based on the same assumptions considered for writing species 
conservation equations for the diffusion layer, the conservation 
equations for species in gas phase and liquid water inside the 
microporous layer can be written as 


-V-: (—Diettm V Cim) = Rw =0 (19) 
-V -Nw + Rw =0 (20) 
where, Di effm is the effective diffusivity and Cim is the concen- 
tration of species i inside the microporous layer. 


3.4. Catalyst layer 


Based on the same assumptions for transport within diffusion 
layer and microporous layer, the conservation equation inside 
the catalyst layer for species 7 in the gas phase can be written as 


=V. (— Di,effr V Cir) + Ri = Rw =0 (21) 


where, Dj effr is the effective diffusivity and Cj, is the concen- 
tration of species i inside the catalyst layer. In Eq. (21), Ri 
represents the consumption of oxygen or production of water 
per unit volume of catalyst layer. Similar to the diffusion layer 
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and microporous layer, Ry in the above equation is applicable 
only for the water vapor balance. Butler—Volmer equation is used 
to define the rate of oxygen reduction reaction 


dal. 
Ro, = = (22a) 
e 
in C F 
dalo = exp( anr ) (22b) 
neF CO, RT eet 


where Co, mem is the dissolved concentration of oxygen in the 
ionomer adjacent to the catalyst site; Co, is the saturation con- 
centration of oxygen inside ionomer. The local overpotential 
(nr) appearing in the Butler—Volmer kinetics is defined by the 
following equation: 


nr = Ad(s, r) — Adge(s, r) 


= {ds E Pr} z {de,s = Per} 


(23a) 
(23b) 


where ġs is the potential of the solid phase and ¢;, is the ionomer 
potential adjacent to active catalyst site. The subscript e denotes 
equilibrium conditions. If standard hydrogen electrode (SHE) is 
treated as the reference electrode, then the solid phase potential 
ds = Veat and Adge(s, r) = Voc, open circuit potential. Therefore, 
the local overpotential (nr) and the ionomer phase potential (¢,) 
are related by the following equation 


Nr = Veat — br — Voc (24) 


Now the question is: how does one determine the dissolved con- 
centration of oxygen inside the ionomer at the active catalyst 
sites? To answer this question, the ideas of macrohomogeneous 
and spherical flooded agglomerate characterization for model- 
ing the catalyst layer region are briefly described in the next two 
sections. 


3.4.1. Macrohomogeneous characterization 

Typically, catalyst layers in PEM fuel cells are prepared by 
synthesizing a uniform mixture of carbon coated with platinum 
nanoparticles and ionomer that is dissolved in a solvent. This uni- 
form mixture is called the catalyst ink. A thin layer of catalyst 
ink is coated on to the membrane and dried. The solvent evapo- 
rates during drying, leaving a thin layer of catalyst coating that 
is usually in the range of 10-50 um thick. Macrohomogeneous 
characterization assumes that any control volume inside the cat- 
alyst layer consists of a uniform mixture of platinum supported 
carbon and ionomer and voids. This is schematically represented 
in Fig. 5. Models based on macrohomogeneous structure assume 
that the concentration of oxygen in the ionomer phase is uniform 
throughout the mixture. Hence, Co, mem in Eq. (22b) is same as 
the concentration of oxygen inside the ionomer at the surface 
i.e. 


Coz,mem = =~ C03, (25) 


where Co, is the concentration of oxygen in the voids inside 
the control volume; Ho, mem is Henry’s constant for oxygen 
between membrane and air. This assumption completely ignores 
the fact that the oxygen concentration will decrease from the 
surface to the active catalyst sites that are embedded inside the 
uniform mixture due to diffusion and reaction. 


3.4.2. Spherical agglomerate characterization 

Spherical agglomerate characterization assumes that the cat- 
alyst particles form agglomerates that are spherical in shape. 
Fig. 6 a shows the schematic of a single spherical agglomer- 
ate in isolation. A thin film of ionomer is also assumed on 
top of the agglomerate [8,7,9]. Another advantage with this 
characterization is that the water produced at the reaction sites 
can be conveniently accommodated using this model. A simple 
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Fig. 5. Schematic of catalyst layer characterized using macrohomogeneous structure. 
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Fig. 6. Spherical flooded-agglomerate for PEM fuel cells 


approach is to assume that the water produced at the reaction 
sites diffuse to the surface of the agglomerate forming a thin 
film before participating in the evaporation/condensation pro- 
cess. A schematic of the spherical agglomerate with a thin film 
of water is shown in Fig. 6 b. The thickness of the water layer 
can be related to the local accumulation and could vary across 
the reaction layer. 

Using spherical agglomerate characterization, one way to 
determine the volumetric rate of oxygen consumption inside 
the catalyst layer is to write species balance equations within 
a spherical agglomerate and solve for Co, mem as a function of 
agglomerate radius. Subsequently, the volumetric rate of oxygen 
consumption can be calculated based on the flux from the profile 
of oxygen concentration inside the ionomer or spherical agglom- 
erate and the number density of spherical agglomerates within 
the catalyst layer. However, this method is very computationally 
intensive [67], as it requires solving for species balance equa- 
tions within a spherical agglomerate for every control volume 
inside the catalyst layer. In order to reduce the time for compu- 
tations, the method of Thiele modulus and effectiveness factor 
was adopted. Moreover, we have also assumed a linear profile for 
the concentration of dissolved oxygen inside the ionomer film. 
Even though the assumption of linear profile introduces some 
inaccuracy for calculating the flux, this particular method has 
been adopted earlier [53,68] and has been shown to be a good 
approximation for predictions. Lin et al. [5] characterized the 
catalyst layer through a cylindrical geometry and have used the 
Thiele modulus and effectiveness factor approach to calculate 
rate of oxygen reduction reaction. Sun et al. [68] have used the 
same approach for the catalyst layer characterized using spher- 
ical agglomerates with a thin film of ionomer. However, they 
have not considered any liquid water effects. The consumption 
of oxygen inside the catalyst layer can be written as 


Ro, = —ġkrxn Co; Ins (26) 
where krxn is the rate constant, given by 


aalo 1 an F 


krxn = 27 
ma = PCR, PORT? (27) 


and ¢ is defined as the effectiveness factor and for spherical 
geometry is given by [66] 


3 
= ee coth(y) — 1) (28) 
The Thiele modulus (W) in the above equation is given by 
k 
v= Fagg A (29) 
O2,mem 


where Fagg is the radius of agglomerate; bg mem 1S the effective 


diffusivity of oxygen in ionomer inside the spherical agglom- 
erates. Concentration of oxygen appearing in Eq. (26) is at the 
ionomer and the spherical agglomerate interface and is given by 


RT cent / HOy,memC0>,r 
{ 1 + (mem/D0,,mem)(1/a1 )okrxn } 


where mem is the thickness of ionomer film; Do, mem is dif- 
fusivity of oxygen in ionomer film; a; is the surface area of 
agglomerate per unit volume of catalyst layer and is given by 
= 2(27 — 2B)(Tagg + mem)” Nagg 
((4/3)7 (Tagg F mem)? Nagg)/(1 — &) 
3 
coe 


E (Fagg + dmem) TT 


(30) 


Co, Ins = 


Me) (31) 


Here 28 is the angle covered by particles or membrane which 
is required for the flow of electrons and protons into the spheri- 
cal agglomerate. The ratio 6/7 is the fraction of surface area 
unavailable for transport of gaseous components and liquid 
water. In the work of [68], the authors have considered only an 
electrolyte film covering the spherical agglomerate. They men- 
tion that a fraction of the agglomerate surface is occupied by 
electrolyte, but they have not tried to quantify the area required 
as pathways for proton transport. In the cylindrical agglomer- 
ate characterization of [53], a liquid water layer on the top of 
the nafion film has been considered. The authors have not con- 
sidered the top and the bottom area of the cylinder for oxygen 
and liquid water transport. An experimental method might be 
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designed to determine such a ratio for a particular system. It 
may also be used as a fitting parameter in a model in absence of 
any experimental data. In order to reduce the number of fitting 
parameters in the model, £ has been treated as zero in this study. 
A sensitivity study will be done in future to study the effect of 
this fraction on the model output. When liquid water forms a thin 
film on top of the spherical agglomerates, it provides additional 
resistance to the transport of oxygen. The additional resistance 
due to water can be accommodated into Eq. (30) by following 
the above analysis, which results in 


(RT ce/ Ho,mem)C03,r 
1+ ((Smem/@1 Doz,mem) 
+ (6w/a1 Do,,w)(Ho,,w/ Ho,,mem))¢krxn 
where ôw is thickness of the water layer; Do,,w is diffusivity of 
oxygen in liquid water; Ho, w is Henry’s constant for oxygen 


between liquid water and air. The thickness of water layer is 
given by 


Coy Ins = (32) 


ErSr 
re E EATE 


Ôw (33) 


ai 
From the stoichiometry of the overall reaction, amount of water 
produced in the cathode catalyst layer is twice the amount of oxy- 
gen consumed. Hence, conservation equation for liquid water 
can be written as 


-V -Nw — 2Ro, + Ry = 0 (34) 


Hydrogen ions produced at the anode catalyst layer travel 
through the membrane and reach the active catalyst sites in the 
cathode catalyst layer. Writing the charge conservation equation 
over a small volume element inside the catalyst layer leads to, 


—V-j,+neFRo, = 0 (35) 


where j, is the local current density inside the cathode catalyst 
layer. For describing the transport of hydrogen ions inside the 
ionomer network, Ohm’s law with an effective ionomer con- 
ductivity and ionomer potential gradient is considered. Hence, 
local current density for the same control volume element can 
be written as 


j: = —Keff V br (36) 
where Keff is the effective conductivity of ionomer inside the 
catalyst layer. It is related to the fraction of volume occupied by 
ionomer (€mem) inside catalyst layer and ionomer conductivity 
(Kmem) by the following equation [5]: 

Keff = ee kinem (37) 
Substituting the Ohm’s law in Eq. (35) leads to 
€nemkmemV “gr + FRo, = 0 (38) 


The cell current density (ice) is calculated by integrating the 
hydrogen ion flux at the cathode catalyst layer and membrane 
interface. It can expressed as 


: | = sel? k Br weel d (39) 
wenl a mem mem ax x=C cell GY 


icell = 


3.4.3. Expressions for €;, €mem and mem 

As mentioned in the introduction, the spherical agglomerate 
characterization is used in deriving the expressions for €mem, €r 
and dmem as a function of the physical and experimental param- 
eters of the catalyst layer such as the thickness (frL) and area 
(arL) of the catalyst layer, weight fraction of platinum on car- 
bon (fpr), catalyst loading (mp), and the weight fraction of the 
ionomer (fmem). The fraction of volume occupied by the voids 
inside the catalyst layer is defined by the following equation: 


Vy 


E&r = (40a) 


URL 


Us 
= 1.0 — — 
URL 


(40b) 


The volume of the solids inside the catalyst layer can be defined 
as the total sum of the volumes of carbon, platinum and ionomer. 
Hence, 


Us = Uc + Vpt F Umem (41) 


Writing the volumes of carbon, platinum and ionomer in terms 
of their mass and density give us 


W 
v=—+F4+— (42) 


Pmem 
The following definitions are used for the weight fractions of 
platinum and ionomer 


Wot 


Spt = (43a) 


Wpt + We 
Wmem 


(43b) 
Wpt + We + Wmem 


fimem = 


Solving Eqs. (43a) and (43b) for we and Wmem, respectively and 
substituting in Eq. (42) leads to 


Wpt [£ 4 1— Sot 
fot | Pot Pc (1 — fmem) mem 


The volume of the catalyst layer (vp) in Eq. (40b) is the product 
of its cross-sectional area (apy) and thickness (tp). Substituting 
for vs and vp in Eq. (40b) gives 


Vs = 


Sie {+s i \ we 
Sot Ppt Pc (1 — fmem)Omem J @RLARL 
(45a) 
1 1— 
S10 f£ Me E ) Ze (45b) 
Sot | Ppt Pc (1 — fmem)Pmem J fRL 


The fraction of volume occupied by ionomer inside the catalyst 
layer is given by 


Emem = e (46a) 
1 
So (46b) 


URL Pmem 


R.M. Rao et al. / Journal of Power Sources 173 (2007) 375-393 


Solving for Wmem using Eqs. (43a) and (43b) and substituting in 
Eq. (46b) leads to 


bee { fmem \ al (47a) 
ee URLPmem | 1.0 — fem) fi pt 
— 1 { fmem \ Mpt (47b) 
tRLPmem | 1.0 — fmem Sot 


Based on a mass balance on Nafion, it can be shown that 


Tage (fmem/1— fmem)(1—(1/mem)) 


3 mem (fpt/ Ppt + — fot/ Pc) + fmem( — fmem)Pmem) 
(48) 


dmem= 


3.5. Boundary conditions, other relationships and 
dimensionality 


The constitutive relationships and other equations defining 
some of the parameters appearing in the above equations are 
given in Table 1. The boundary conditions for different sections 
of the cathode are given in Table 2. Gas entering at the cathode 
inlet is assumed to be a mixture of air and water vapor with 
no liquid water saturation (BC1). At the boundaries of the gas 
flow channel/Diffusion layer, diffusion layer/microporous layer, 
microporous layer/catalyst layer, continuity of the variables and 
flux continuity are imposed. Boundary conditions 2 and 3 sug- 
gest that the liquid water saturation is continuous at the gas flow 
channel/diffusion layer and the liquid water flux is also equal. 

The boundary conditions for the liquid water saturation (s) 
need a detailed discussion in this context. The conservation equa- 
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tion for liquid water in gas flow channel is first order in y. We have 
considered the inlet saturation of gas as the boundary condition 
for this equation. The conservation equations for liquid water 
in the gas diffusion layer (GDL), microporous layer (MPL) and 
catalyst layer (CL) are second order in x. The obvious and well 
accepted boundary conditions are continuity of flux at CL/MPL 
interface and MPL/GDL interface. No flux condition is assumed 
at CL/membrane interface. But three more boundary conditions 
are required—each at channel/GDL interface, GDL/MPL inter- 
face and at MPL/CL interface. 

A survey of the published literature reflect the different 
boundary conditions adopted by different authors. Pasaogullari 
and Wang [4] considers the catalyst layer to be ultrathin and 
therefore the oxygen reduction reaction is assumed to take place 
at the interface of PEM and diffusion layer. Therefore, only 
the continuity of flux is used as a boundary condition at this 
interface. In their work, they assume the liquid saturation at the 
channel/GDL interface to be zero. In another work by [44], the 
same boundary conditions as above have been considered. In the 
work of [24], the flux of liquid water at channel/GDL interface 
has been considered zero at the inlet. At the outlet, the authors 
have considered ds/dy = 0 as a boundary condition. The authors 
have also considered the catalyst layer to be ultrathin. Therefore 
the single boundary condition of continuity of flux is sufficient 
for them at GDL/membrane interface. Natarajan and Nguyen 
[69] assume that once the gas stream in the gas flow channel is 
saturated, the boundary conditions of liquid water saturation for 
subsequent volume elements at channel/GDL interface become 
0.1. In that work, they have considered reaction layer as ultrathin 
requiring only the flux boundary condition. In the work of [53], 


Table 1 
Constitutive relations 
Parameter Expression Reference 
Di eff Gu al — s)?/ 2 Dim Bruggeman relation 
Ky(s) Kwos [53] 
Mya 
da P 
ÍRL 
1.0 
Cc 
02 Ho, mem 
666.0 
Ho, mem 1.33 exp | — 11] 
cell 
498.0 
Ho,,w 5.08 exp (- ) 72] 
~ cell 
2768 
Do, mem 3.1 x 1077 exp (- ) Femti 
4/9 Teell 
aa el? Do, ,mem Bruggeman relation 
1 1 
100.0(0.005139A,, — 0.00326 1268.0 | ——— — 36 
Kmem ( w ) exp | (sa Teel )| ] 
day 0.043 + 17.8lay 39.8542, H 36.043, for0 < aw < 1, 14.0 + 1.4(ay — 1) forl < aw < 3, 16.8 for aw > 3 36] 
Pwr 
aw pst 
WwW 
1.02 
sat T —4.9283 1023-5518—(2937.8/T)) 
Pw 1000 
1.98 + 0.03242, ) 3 
—— 10° 74 
Pmem ( I+ 0.06482, [74] 
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Table 2 
Boundary conditions for PEM fuel cell cathode shown in Fig. 2 


Number Location 


Boundary condition 


Comments 


Gas flow channel 
BCI Y=0 


Diffusion layer 


Ci = Cio, Sge = 0 


i = O2, No, H20, liquid saturation at 
inlet 


BC2 X=0,VY Cia = Ci gc, Sge = Kisa Continuity at surface 

BC3 X=A,VY Di eiim 2 = — Dieta it, Nwmle = Nwalx Flux continuity in X direction 
BC4 0<X<A,Y=0 Di effa wis = 0, Nwaly = 9 Zero flux condition in Y direction 
BCS 0<X<A,Y=L Di effa ti = 0, Nwaly = 0 Zero flux condition in Y direction 


Microporous layer 


BC6 X=A,VY Cia = Cim, Sa = K25m Continuity at surface 
BCT X = B,VY Diettm SE = — Diet GE, Nwmle = Nwrlx Flux continuity in X direction 
BC8 A<X<B,Y=0 — Di effîm n = 0, Nwmly = 0 Zero flux condition in Y direction 
BC9 A<X<B,Y=L — Di effîm ia = 0, Nwmly = 0 Zero flux condition in Y direction 
Catalyst layer 
BC10 X = B, YY Cir = Cim, Sm = K35r, E =0 Continuity at surface in X direction, 
zero flux for Ht ions 
BC11 X=C,VY Di effr Cis = 0, Nwrlx = 0, 6 = 0 Zero flux condition in X direction, 
reference point 
aCix 3/2 abr se on oe 
BC12 B<X<C,Y=0 Dietir >> = 9, Nwrly = 0, —EmemKmem py = 0 Zero flux condition in Y direction, 
f i zero flux for Ht ions 
BC13 B<X<C,Y=L Diettr ois = 0, Nwrly = 0, er enate =0 Zero flux condition in Y direction, 


zero flux for HT ions 


the transport of liquid water in the gas flow channel is not consid- 
ered in their one-dimensional model. The liquid water saturation 
is considered to be same at the inlet and along the channel. In 
this model, liquid water in the reaction layer has been modeled. 
The liquid water conservation equation in the GDL and CL is 
second order in s. But only one boundary condition is mentioned 
in their work. 

In the subsequent work of [70], the same boundary condi- 
tions as before have been considered. The boundary condition 
that may be used at these various interfaces are the relation 
between interface saturations. But the interface saturation is a 
very complicated function of the surface properties and the oper- 
ating conditions. One of the most comprehensive discussions 
about the interface saturation is first given by [30]. The authors 
expressed the interfacial saturation at channel/GDL interface 
as a function of channel gas velocity, contact angle and local 
current density. Although liquid water generation is related to 
the local current density, interfacial transfer of liquid water to 
vapor phase has been considered in the present model. Therefore 
a part of the liquid water gets transferred to vapor phase as it 
flows across the reaction layer, microporous layer and the diffu- 
sion layer to the flow channel. As long as the saturation pressure 
of water at the operating condition is not exceeded, the liquid 
water that has been transferred to vapor phase remains in vapor 
phase. Therefore it is unlikely that it should affect the saturation 


at the interface. So in the view of the authors: 


at X=0,VYY Sge = fi(Sa, 9, uy) (49) 


A detailed experimental test is required to determine such a 
function for a particular system. A detailed procedure can be 
found in the work of [71]. In this model, for simplicity , a linear 
relationship between Sgc and Sq is assumed 


at X=0,VY Sge = Ki Sa (50) 


where Kı is a constant. The constant is always less than one 
and might be used in either side of the equation to maintain the 
saturation s (which is a volume fraction) to be less than one. 
Similarly the saturation at the GDL/MPL interface will depend 
on the wettability of both the surfaces. So 


at X = A, VY Sos = fo(Sm, 01, 02) (51) 
where 6 and 62 are the contact angles of backup substrate and 
the diffusion layer respectively. If the medium is hydrophilic, 
0 < 90° whereas for a hydrophobic medium, @ > 90°. The 
velocity of gas in the porous media may also affect the inter- 
face saturation. In this work, the gas phase velocity has been 
considered negligible. In the absence of any experimental data, 
and for simplicity, here also a linear relationship between Sg and 
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Sm is assumed: 
at X = A, VY Sbs = KoSm (52) 


Ka is a constant in Eq. (52). Because of similar reasons, one of 
the boundary conditions at the MPL/CL interface has been taken 
as: 


at X = B, YY Sa = K3S, (53) 


K3 is a constant here. In this context, it may be mentioned that 
as (dP, /ds) has been assumed constant in this study, the depen- 
dence of P, on liquid saturation at different contact angles has 
not been explicitly modeled. The dependence has rather been 
assumed constant in this operating range and for the given sur- 
face wettabilities. Although the model does not consider the 
effect of varying hydrophobicity on the interphase saturation, 
the use of different constant values for different layers improves 
the predictive capability of the model significantly. 

Fig. 2 depicts the geometry of the gas supply by using a 
machined graphite plate with straight flow channels. As we can 
see that the channels are separated by ribs or shoulders, which 
occlude half the flow area on top of the backup substrate. This 
necessitates considering flux variation in all the three directions 
(X, Y and Z) inside the porous regions. The ribbed geometry of 
the graphite plate causes the local flux to increase and causes 
a greater concentration drop than for a uniform flow. The idea 
of considering a uniform flow with a thicker diffusion layer that 
drops the concentration to an equivalent amount was expressed 
by [75]. They compared different geometries with various flux 
and concentration profiles for both steady state and dynamic 
cases and suggested an increase in the thickness of the diffusion 
layer by a factor of 0.6. In order to reduce the dimensionality 
in the steady and the subsequent dynamic models, the thick- 
ness of the diffusion layer has been increased by a factor of 
0.6. Hence, in the steady state and the dynamic models, spa- 
tial variations are considered only in the Y direction inside the 
gas flow channel and in X and Y directions inside the porous 
regions. 


4. Modeling methodology 


The steady state model described in the previous section can 
be solved in two different ways. One method is to solve the 
model equations by considering cathode voltage, design and 
other model parameters as inputs and predict the cell current or 
current density. The second method is to treat the cell current 
or current density along with the design and model parameters 
as inputs and predict the cathode voltage. If we study the model 
equations carefully, we can see that the voltage appears inside an 
exponential term. Hence, it may be computationally less expen- 
sive to solve the model using the first method, i.e., treat voltage 
as an input and predict the current density. Preliminary trials in 
solving the model equations using both methods confirmed this 
point. 

The model equations are essentially partial differential equa- 
tions that are non-linear due to the interactions between the 
different variables. Since it may not be possible to obtain an 


analytical solution, the system of equations have been solved 
numerically. The next section describes in detail the methodol- 
ogy that was adopted for modeling and solving the system of 
equations. 


4.1. Model development 


The model development was carried out using Maple and 
subsequently, MATLAB was used to solve and simulate the 
model. The partial differential equations (PDEs) and the bound- 
ary conditions corresponding to the steady state model are setup 
in Maple in a symbolic manner. They are discretized in spa- 
tial variables using finite difference techniques and converted 
into a system of non-linear algebraic equations (NAEs). The 
discretization was done using an in-built Maple library func- 
tion called convert and a Maple procedure, which can be used 
to produce finite difference approximations for PDEs. The pro- 
cedure takes an expression and specifications for making finite 
difference approximations such as, backward, forward, or cen- 
tral difference formulae and the desired order as inputs. The 
output from the procedure is an expression representing the 
finite difference formula with appropriate indices. Hence, the 
PDEs and boundary conditions are converted to a list of NAEs 
in Maple. Further, the list of NAEs is converted into a proce- 
dure in Maple using an in-built library function called proc. 
Using a code generation package in Maple, the Maple proce- 
dure is converted into a MATLAB function, which essentially 
consists of all the equations. Finally, the MATLAB function can 
be called by an equation solver by passing appropriate parame- 
ters. In the present case, the non-linear algebraic equations are 
solved using fsolve, which is a part of MATLAB s optimization 
toolbox. The whole procedure is schematically represented in 
Fig. 7. 


4.2. Computational issues 


Typically, there are several computational issues that are 
encountered while solving systems of PDEs and NAEs. Some 
of the issues and the techniques that have been used and/or 
developed to overcome them are highlighted below: 


e Scaling of variables: several trials in solving the steady state 
and dynamic model equations demonstrated that the scaling 
of the variables is important due to the existence of highly 
non-linear terms in the equations. In its unscaled form, the 
system of equations could not be solved for the whole range 
of operating voltage. Hence, the variables in the equations 
have been appropriately scaled using the inlet conditions as 
reference values. 

e Calculation of Jacobian: During preliminary trials in solv- 
ing the model equations, it was observed that the numerical 
calculation of Jacobian was computationally very expensive. 
To overcome this, analytical expressions of non-zero ele- 
ments in the Jacobian matrix were generated in Maple. These 
expressions were subsequently used by the solver fsolve while 
solving the system of equations. The use of analytical Jaco- 
bian greatly reduced the computational effort. 
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Fig. 7. Flowchart for modeling methodology. 


e Grid sensitivity: Another important issue that has been studied 
is the sensitivity of the model predictions to the number of 
grids considered for writing the finite difference formulae. 
Different number of grids were considered in the Y direction in 
gas flow channel (GC), both the X and Y directions in diffusion 
layer, the microporous layer and the catalyst layer (RL). The 
steady state model was simulated and the sensitivity of the 
i-v curves is illustrated in Fig. 8. Based on this an appropriate 
number of grid points were chosen. 


Mesh: GC 11; BS 11 x 11; DL 11 x 11; RL 11 x 11 
Mesh: GC 16; BS 11 x 16; DL 11 x 16; RL 11 x 16 
Mesh: GC 16; BS 16 x 16; DL 16 x 16; RL 16 x 16 
Mesh: GC 11; BS 11 x 11; DL 11 x 11; RL 21 x 11 
Mesh: GC 16; BS 11 x 16; DL 11 x 16; RL 21 x 16 
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Fig. 8. Steady state polarization curves vs. mesh density. 


5. Steady state model validation 


Experimental data (i-v curves) for a PEM fuel cell was used 
to validate the steady state model. It was shown by [76] that 
experimental validation of multiphysics PEFC models must be 
done against data at the distribution level. However, for the scope 
of this study, the steady state model is validated using the overall 
i-v curves. Further validation of these models with distributed 
data is currently under progress. 

The geometry of the PEM fuel cell is similar to that illus- 
trated in Fig. 2. The cell consists of eight straight channels in 
the graphite plate and an effective MEA area of about 20cm’. 
The cell is operated at 65 °C. The design and other input param- 
eters associated with the experimental data are given in Table 3. 
The fuel cell was operated with three different air flow rates. 
For every flow rate, the cell characteristic data was generated by 
operating the cell at different currents in the range of 0-10 A. 
Cathode voltage measured against standard hydrogen electrode 
(SHE) was recorded for each operating current. The steady state 
polarization data for the three different air flow rates are shown in 
Fig. 9. The data is an average of three different runs for every one 
of the flow rates. The next section describes about how the model 
described in Section 3 was validated against the experimental 
data represented in the above figure. 


5.1. Model validation 


As one can see from the description given in Section 3, 
various assumptions were introduced for writing the model 
equations. The validity of the assumptions are verified by com- 
paring the corresponding simulation with the experimental data. 
The following modeling assumptions were tested: (i) transport 


Table 3 
Inputs and design parameters for experimental data 
Parameter Base case value Units Comments/reference 
Inputs 
Poat zx 1.0 atm 
Pag 0.2, 1.5 and 2.0 Ipm 
Veat 0.95 to 0.25 Vv 
Teel 338.15 K 
Tair 343.15 K 
RHair 100.0 % 
Design parameters 
Weh 0.002 m 
hen 0.002 m 
L 0.07 m 
GDL 350.0/0.6 x 1076 m 
MPL 450.0 x 1076 m 
IRL 45.0 x 1076 m 
Nch 8 
Weell 0.03 m 
Mpt 0.2 mg Ptcm~?RL7! 
apt 20 m? Ptg! Pew! 
Jot 0.20 
fmem 0.25 
Ebs 0.70 
€d 0.50 
€r 0.50 
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Fig. 9. Steady state polarization data for PEMFC cathode. 


governed by diffusion in the bulk inside the porous regions 
(diffusion layer, microporous layer and catalyst layer), (ii) trans- 
port governed by simultaneous bulk and Knudsen diffusion, 
(iii) cathode catalyst layer characterized using a macroho- 
mogeneous structure, (iv) cathode catalyst layer characterized 
using a spherical agglomerate structure, with and without liq- 
uid water effects. Various models, each representing the above 
assumptions were developed using the methodology described 
in Section 4. For model validation, the following procedure was 
adopted. The steady state experimental i—v data corresponding 
to Fair = 1.5 lpm was used to compare and tune the model pre- 
dictions. The base case model parameters used for predicting 
the steady state polarization curve are listed in Table 4. 

Fig. 10 shows the model predictions against the experimental 
data for the case where the cathode inlet air flow rate was 1.5 lpm. 
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Fig. 10. Model predictions vs. experimental data. Fajr = 1.5 lpm. 


First, a macrohomogeneous model was used with the assump- 
tion that the transport in the porous regions take place only due 
to diffusion in bulk. The predicted i—v is shown using ‘square’ 
legend. It can be observed that the model prediction matches 
experimental data only in a small range of voltage, 0.95-0.75 V. 
Below 0.75 V, model predictions of current density are much 
higher than the experimental values. On the other hand, con- 
sidering Kundsen and bulk diffusion in parallel, current density 
values predicted by the model (shown in open circles) are much 
less than the experimental values. Model studies where transport 
in the porous regions was assumed to be purely due to Kundsen 
diffusion were also studied and the predicted current density 
values were found be lower than experimental values. At the 
end of this exercise, it was concluded that transport in porous 
regions must be predominantly occurring due to diffusion in 


Table 4 
Model parameters for the base case 
Parameter Base case value (final value) Units Comments/reference 
Model constants 
F 96485.3 Cg7leq. 
Ne 4 
R 8.314 Jg-! mol! K7! 
Pe 1.8 x 10° kg m~? 
Ppt 21.45 x 103 kg m~? 
Pw 983.21 kg m3 
Model parameters 
Tage 1x 1077 m Middleman [7] 
ig 3.72 x 1075 (3.84 x 1073) if Vee > 0.79 V Am? Pt! Parthasarathy et al. [77] 
1.1 x 107? (3.2 x 107!) if Vee < 0.79 V Am? Pt"! Tuning parameter 
a 1.0 if Veen > 0.79 V Parthasarathy et al. [77] 
0.625 if Veen < 0.79 V 
ke 100.0 s7! For all layers [53] 
ky 100.0 atm™! s7! For all layers [53] 
Kwo 7x 107! m? Catalyst layer [53] 
7x 10-3 m? Microporous layer [53] 
7 x 1078 (1 x 107!) m? Diffusion layer (tuning parameter) 
Ki 1.0 
Ka 1.0 
K3 1.0 
(—dP,/ds) 113.68 N m~? Catalyst layer [53] 
56.84 N m~? Microporous layer [53] 
56.84 (28.42) N m~? Diffusion layer (tuning parameter) 
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the bulk. Hence, all subsequent models considered only bulk 
diffusion. 

Next, a cathode model in which the catalyst layer was char- 
acterized using spherical agglomerate structure was simulated 
to predict the i-v data. The radius of the agglomerate consid- 
ered for the model is shown in Table 4. In this model, it was 
assumed that all the water that is produced inside the cathode 
catalyst layer is in vapor form. The corresponding i—v data is 
shown using “x” legend in Fig. 10. It can be observed that 
even though the predicted current density values are higher than 
experimental data, they are significantly lower than the corre- 
sponding values predicted by the macrohomogeneous model. 
Subsequently, the spherical agglomerate model was extended 
to include the effects of liquid water in the porous regions, as 
described in Section 3. With liquid water effects, the model pre- 
dictions matched very well with experimental data, as shown by 
the curve using “*” as legend. The match between experimental 
data and model prediction is clearly shown in Fig. 11. 

Hence, it can be seen that when the catalyst layer is character- 
ized with spherical agglomerates and transport in porous regions 
is by diffusion in bulk, the model predictions match very well 
with experimental data. The macrohomogeneous model assumes 
that within a small control volume in the catalyst layer, the dis- 
solved concentration of oxygen is uniform throughout the homo- 
geneous mixture of ionomer and platinum supported carbon. 
However, in the model based on spherical agglomerates, the dis- 
solved concentration decreases from the surface of the ionomer 
to the active catalyst sites. Table 4 shows some of the model 
parameters that were used as tuning parameters. The final values 
of the tuning parameters are shown in parentheses in the table. 

To check the validity of the model, experimental data corre- 
sponding to two very different air flow rates of 0.2 Ipm (low) 
and 2.0lpm (high) were compared to the model predictions. 
The results are shown in Figs. 12 and 13, respectively. It can 
be seen that the model predictions are close to the experimen- 
tal data. Hence, with the given set of experimental data it can 
be concluded that the cathode model predictions are valid to a 
large extent with small deviations in some range of the operating 
voltage. Some of the assumptions that were made in Section 3 
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Fig. 11. Spherical agglomerate model predictions vs. experimental data. Fair = 
1.5 Ipm. 
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Fig. 12. Model predictions vs. experiment. Fair = 0.2 lpm. 


for transport inside the gas flow channel and the porous regions 
were motivated by the cell design for which the experimen- 
tal data was collected. The proposed cathode model has good 
predictive capabilities for small scale PEM fuel cells. The pre- 
dictions of the proposed model needs to be tested against larger 
fuel cell designs. 

Depending on the wettabilities of the materials at the different 
layers and operating conditions of the cell, the liquid water satu- 
ration at the interphases may vary widely. Meng and Wang [30] 
performed a sensitivity study on the value of the parameter that 
they use in their study at the channel/GDL interface. The authors 
have considered two values of the parameter at the channel/GDL 
interphase and have observed how a decrease in the interfacial 
liquid saturation improves the cell performance. In the present 
work, we have observed the effect of the interfacial saturation at 
MPL/CL interface on the cell performance. Fig. 14 shows that 
as K3 varies from 1.0 to 0.4, the performance of the cell does 
not change much. As the value of K3 decreases further, the cell 
reaches its limiting current density at a higher cell voltage than 
before. At a value of 0.2, the cell current density gets limited to 
about 200 mA/cm?. As the value of K3 gets reduced, the out- 
flow of the liquid water from the cell gets reduced because of 
the reduced driving force. In other words, as the saturation level 
in the catalyst layer increases, the thickness of the water film 
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Fig. 13. Model predictions vs. experiment. Fair = 2.0 lpm. 
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Fig. 14. Polarization curves for different values of K3. Fair = 2.0 lpm. 


around the spherical agglomerate increases, pushing the cell to 
a mass transfer limited regime. A similar study is carried out 
at the low flow rate of 0.2 lpm also. The nature of the plots do 
not change, but the cell reaches its limiting current density at a 
higher cell voltage than that of the simulations at 2.0 lpm for any 
given value of K3. The simulation results clearly show that the 
interfacial liquid coverage plays a very important role in deter- 
mining the limiting current density of a cell in a given operating 
condition. 


5.2. Utility of the steady state model 


The steady state cathode model can be used to study the effect 
of various inputs and design parameters (Table 3) on i-v char- 
acteristic curves. Inputs such as, flow rate of air, inlet pressure, 
inlet air humidity, inlet air temperature, and operating cell tem- 
perature can be considered. The effect of inlet air flow rate on 
i-v characteristic curve was already seen in Figs. 10, 12, and 13 
during model validation. In addition, effect of design param- 
eters such as, dimensions of the cell, catalyst loading (mp), 
weight fraction of platinum on carbon (fp), and weight frac- 
tion of ionomer inside catalyst layer (fmem) can also be studied. 
Moreover, the steady state model can be used to study the spa- 
tial distribution of partial pressure or concentration of species, 
liquid water saturation, and ionomer phase potential inside dif- 
ferent regions of the cathode. The effect of the operating cathode 
voltage on spatial distributions inside the catalyst layer is shown 
in Figs. 15-17 . Fig. 15 shows the spatial distribution of partial 
pressure of oxygen inside the catalyst layer for three different 
cathode voltages: Veat = 0.70, 0.50, 0.30 V when inlet air flow 
rate is maintained at 0.2 lpm. For X direction inside catalyst layer, 
number “1” corresponds to microporous layer and the catalyst 
layer interface (X = B in Fig. 2), and number “11” corresponds 
to the catalyst layer and the membrane interface (X = C in 
Fig. 2). For Y direction inside the catalyst layer, number “1” 
corresponds to the inlet side (Y = 0 in Fig. 2) and number “11” 
corresponds to the outlet side (Y = L in Fig. 2). 

Fig. 16 shows the liquid water saturation inside the catalyst 
layer for the same three cathode voltages. We can see that when 
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Fig. 15. Spatial variation of partial pressure of oxygen inside catalyst layer. 
Fair = 0.20 Ipm. 
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Fig. 16. Spatial variation of liquid water saturation inside catalyst layer. Fair = 
0.20 Ipm. 


the cathode voltage is low (0.3 V), implying high potential driv- 
ing force for the chemical reaction, there are regions inside the 
catalyst layer where 70% of the voids are filled with liquid water. 
The region inside the catalyst layer (number “1” along Y direc- 
tion and number “11” along X direction) corresponds to control 
volumes near the membrane and the inlet. Fig. 17 shows the 
spatial variation of ionomer phase potential inside the catalyst 
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Fig. 17. Spatial variation of ionomer potential inside catalyst layer. Fair = 
0.20 Ipm. 
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layer. One can see that the gradients along the X direction become 
steeper when cathode voltage decreases from 0.70 to 0.30 V. It 
can be observed that all the three surface have a common edge at 
the top. This is due to the fact that the catalyst layer is interfaced 
with a standard hydrogen electrode for voltage measurements 
and this edge is taken as the reference. 

Steady state models can be used to do optimization studies 
wherein some inputs or design parameters are allowed to vary 
within a range of values to maximize or minimize an objective 
and at the same time satisfying one or more criteria. In a related 
work [10], the steady state cathode model was used to study an 
optimization problem where the objective was to maximize the 
cell current density. In order to do this, catalyst loading (mp), 
weight fraction of platinum on carbon (fpt), weight fraction of 
ionomer inside catalyst layer (fmem), and catalyst layer thickness 
(tpL) are treated as optimization variables. The details of the 
optimization formulation, numerical procedure and results are 
described in [10]. 


6. Conclusions 


A two-dimensional steady state model for the cathode of a 
PEM fuel cell was presented in this work. Various modeling 
assumptions and catalyst layer characterizations were tested dur- 
ing model validation. It was shown that the model predictions 
using the spherical agglomerate characterization of the catalyst 
layer best fits the experimental data. An approach to correlate 
experimental conditions to model parameters was suggested. 
Using this correlation, it is possible to perform multi-parameter 
optimization studies for fuel cell systems. The results of such 
an optimization study is presented in [10]. Improvements to the 
model proposed in this paper can lead to the use of such mod- 
els in design of fuel cell materials for better performance. It was 
also shown that the wettability of the different layers plays a very 
important role in getting the best performance out of a PEM fuel 
cell. 
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